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Abstract 

A new formulation (including the 
choice of variables, their non- 
dimensionalization and the form of the 
artificial viscosity) is proposed for the 
numerical solution of the full Navier- 
Stokes equations for compressible and 
incompressible flows with heat transfer. 

With the present approach, the same 
code can be used for constant as well as 
variable density flows. The changes of the 
density due to pressure and temperature 
variations are identified and it is shown 
that the low Mach number approximation 
is a special case. At zero Mach number, the 
density changes due to the temperature 
variation is accounted for, mainly 
through a body force term in the 
momentum equation. It is also shown that 
the Boussinesq approximation of the 
buoyancy effects in an incompressible 
flow is a special case. 

To demonstrate the new capability, 
three examples are tested. Flows in driven 
cavities with adiabatic and isothermal 
walls are simulated with the same code as 
well as incompressible and supersonic 
flows over a wall with and without a 
groove. Finally, viscous flow simulations of 
an oblique shock reflection from a flat 
plate are shown to be in good agreement 
with the solutions available in literature. 

Introduction 

In a previous work[l], the authors 
proposed a formulation for both 
compressible and incompressible viscous 
flow simulation. First the density is 
eliminated in terms of the pressure and the 
temperature via the perfect gas equation 
of state. This step by itself is not sufficient 
simply because the equation of state is not 
valid for incompressible flows. The 
formulation is completed using the 


perturbation of the pressure and the 
temperature relative to reference values as 
the dependant variables. It is shown that 

the density in terms of these new variables 
approaches a constant as the reference 

Mach number vanishes. The above 
formulation is generalized in the present 

paper to allow for incompressible flows 
which are not necessarily isothermal. 

To obtain a numerical solution, an 
artificial dissipation is introduced by 
adding to the governing equations the 
Laplacians of the pressure and the velocity 
components. An improved model is also 
tested which is based on a partial least 

square procedures. The continuity 
equation is modified by a Poission's 

equation for the pressure similar to that of 
Harlow and Welch[2], and Harlow and 

Amsden[3]. The momentum equations are 
also modified by Poission's equations of the 
velocity components. The first 
modification is obtained by taking the 
divergence of the momentum equations, 
while the second modification can be 
related to a vector identity relating the 
Laplacian of the velocity vector to the 
gradient of its divergence and the curl of 

the vorticity. In both modifications, the 
evaluation of the nonhomogeneous terms 
of the Poissions equations are lagged as in 
the deferred correction procedures. 

The energy equation is augmented 
with second order terms of the total 
enthalpy obtained via minimizing the 
squares of the convective terms. This 
modification is very small in the 
neighborhood of a solid surface and can be 
interpreted as an artificial streamline 
diffusion as in the work of Hughes etal[4]. 

The present numerical solutions are 
obtained using a standard Glarkin 
procedure. The resulting nonlinear system 
of equations are solved via Newton's 
method. At each iteration a direct solver 
based on banded Gaussion elimination is 
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employed. The use of finite element 
discretizations and direct solvers are not 
necessary to obtain numerical solutions 
based on the present formulations, and- 
other viable alternatives are, for example, 

finite volumes and iterative procedures. 

In the following, the derivation of the 
governing equations and the applications 
to some test problems are discussed. 

Governing Equations 

For steady compressible viscous flows, 
the continuity, momentum and energy 

equations can be written in terms of the 

primitive variables (p, p, ~q ) including the 
effects of a body force as : 

V. p^j =0 
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V. pq q = V : x + p f 

V. p'q H = V.kVT + V.(-tiq) + p7.'q (1) 
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For convenience it is assumed that X = - 
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2 y. and p = p R T. The two constants R and y 

are related to the specific heat constants 
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c and c ; R = c - c and y = ~ . For the 
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P 

derivation of the above equations, see, for 
example, Liepman and Roshko [ 5 ] 

A standard non-dimensional form is 
obtained using the reference values of p 
and q in the far field of external flow 
problems. The pressure is usually 
2 

normalized by p q and the temperature 
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by q / c . If L is a characteristic length, 
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two parameters appear in the equations 
namely the Reynolds and Peclet numbers, 
where 


Re =- 


p q L 
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The Peclet number is the product of 
Prandtl and Reynolds numbers, where 


Pr = 


k /c 
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Equations ( 1 ) becomes 
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V. pqH= — V. kVT+— V .(t .q)-— q . pK 

(5) 

In equations ( 5 ), the relative effect of 
gravity is identified by the Froude number, 
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g L 
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In natural convection problems, the 
variation of density due to temperature 
difference AT creates a buoyancy term in 
the momentum equation. To first order 
accuracy, the density variation would be p 
S P ( 1 - p AT) where P is the thermal 


expansion coefficient, hence the buoyancy 
term is given by 


Gr £ 
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where 
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The above formulation is not suitable 
for incompressible flows since in the limit 
of zero free stream Mach number, both the 
normalized pressure and temperature are 
unbounded. Two new variables were 

introduced in the previous study to avoid 

* * 

this problem, namely p and T where 
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T = T - T = T 


(y-1) iv? 

oo 

Hence, the equation of state gives 
y N? p +1 


( 10 ) 


P = 


( 11 ) 


(y-l)N? T +1 
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As the Mach number vanishes, the 

normalized density p approaches 1.0 , i.e. 
the reduced incompressible flow is 
isothermal. To allow for density variations 
due to temperature changes in the 
incompressible limit, in cases of adiabatic 

walls as well as walls with specified 

* 

temperatures, the variable T is replaced 
by 


- T 
T_ T 
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Equation (11) becomes 
2 * 

yM p +1 


P = ' 


( 12 ) 


(13) 


In the limit of zero Mach number, p T 
approaches 1 and the proper general 

dependance of the density on the 
temperature is recovered. The isothermal 

flow is of course a special case of the above 
relation. 

The continuity and the momentum 
equations are unaltered, the energy 

equation becomes 


--■=?“ l ~ Ec- -■* Ec - ^ — * 

y pqH= Pe v kVT+ Re V(T,q) - q P K 
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where 
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H=T+ Ec q /2 
and, 

2 
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Ec = — = (y-1) M (Ec is the Eckert 
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number) 

The incompressible limit of Equation (14) 
is 


V. p^T=^kVT (15) 

Here, the temperature ratio T /T , 

V/ 00 

where T is the average wall specified 
w 

temperature, enters only through the 
boundary conditions. 

Alternatively, one can choose 
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and, in this case 
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Thus, the present formulation is valid 
for compressible and incompressible flows, 
with adiabatic or specified temperature 
walls. Moreover, it is clear from equation 

(13) or (13 ) that the low Mach number 
approximation ( see for example Rhem and 
Baum [6], Majda[7] and Markle[8]) is a 
special case. The Boussinesq approximation 
of the buoyancy effects in an 

incompressible flow is also a special case of 
the present formulation. 

It should be mentioned that the above 
formulation is not restricted to perfect 
gases. A more general equation of state can 
be written in the form 
dp dp 

p =p .*aT aT< '5p ap 


(16) 



or in the non-dimensional form 

__ dp ~ dp * — 

p = l+-fAT+^4p (16) 
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The last term of equation (16) always 
vanishes in the limit of zero Mach number. 

Numerical Method 

It is well known that centered schemes 
permit, in general, odd and even 
decoupling of the discrete pressure 
field[9]. To avoid this problem, different 
interpolations are used for the velocity and 
the pressure in the standard finite element 
analysis of incompressible flows [10], [11], 

Recently, Pironeau [12] addressed this issue 
for compressible flows as well. 

It is also well known that centered 
schemes produce oscillatory solutions of 
convection-diffusion equations with high 
Reynolds numbers, unless impractical 
excessively fine meshes are used 

In the present study, artificial 
dissipation is introduced explicitly in all 
equations, to eliminate the wiggles and to 
allow for capturing shocks and contact 
discontinuities. Two forms of artificial 
viscosity are considered. In the first 
method, the governing equations become 
2 

V p q = Ej V p 
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V pqq+Vp V: t+— P K=e 2 V q 

\fpqH-^kVT-|^V.(T% + pf.g=e 3 V H 

(17) 

where e's are small parameters of the order 
of the mesh size. A standard Galerkin finite 
element method is applied to calculate the 
solution of equations (17). This form has 
been investigated before for both 
incompressible and compressible flows 
with the standard separate formulations. 
With the present unified approach, the 
same code is used to calculate compressible 
and incompressible “cases. 

In the second method, the Laplacian 
terms are balanced with nonlinear terms 
obtained by manipulating the original 


equations. For example, the momentum 
equations can be written in the form 

V p = g (18) 

A poission's equation is constructed by 
taking the divergence of equation (18) and 
allowing a variable (positive) artificial 
viscosity coefficient, one arrives at 

V,eVp=V.eg (18) 

Equation (18 ) can be a lso obtained 
from minimizing the functional j e (V p - 

— > •• . 

g ) , (V p * g ) with respect to p, assuming 
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g is known. The continuity equation is 
then modified by the Poission's equation 

(Tg): 

Similarly, the Laplacian of the velocity 


components can be balanced 

using the 

vector identity 


2_> -> 


V q =VS- VXco 

(19) 


where 


S = V 
a> = V X q 

To allow for a variable viscosity 
coefficients, one can minimize the 

— > 

functional (assuming S and co known) 

J e (V ."q - S ) 2 + e (V xlf -w) 2 
with respect to to obtain 

V ,e Vq = VeS-VX £0 + f(Ve) ( 1 9) 

In equations (18 ) and (19 ), the 
— > — * 

quantities g , S and to are obtained at each 
node from their definitions using a 

standard Galerkin finite element method 
with the same interpolation used for the 
other variables. The evaluation of these 
terms are lagged and their contributions to 
the Jacobians are neglected. 

For the energy equation, the 
modification is obtained via minimizing 

- 2 

the functional i £ (p q V H + it) with respect 
to H. For convenience, it is dropped with the 
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justification that the artificial dissipation 
is mostly needed in the inviscid adiabatic 
part of the flow where the relatively 
coarse mesh is not capable of resolving the 
Jt term. In the neighborhood of a solid 
surface, a fine mesh is required anyway 
and there is no need there of artificial 
viscosity terms. The present modification is 
very small there since it is scaled with the 
velocity. The same remark is applicable for 
the treatment of the momentum equations 
where the viscous stress terms can be 

— ) 

ignored in the evaluation of the g term. 

The variational formulation of the 
artificial dissipation terms provides a 
natural treatment of the numerical 
boundary conditions. Upon integration by 
parts, the resulting line integrals are 
simply ignored. 

Because the modification terms, for the 
continuity, the momentum and the energy 
equations are obtained separately by 
adopting a partial least squares procedure 
for each case the resulting algebraic 
system of equations are not necessarily 
symmetric.. A full least squares procedures 
for all equations coupled together has been 
successfully used by the first author to 
introduce dissipative terms for the solution 
of Euler equations simulating transonic 
flows with sharp shock waves [13]. In this 
case, it is possible however to construct a 
symmetric positive definite system at each 
Newton’s iteration by a proper choice of 
the coefficients of the artificial terms [ 
Jiang& Povinelli[14], Unfortunately, such 
choices result in smearing the 

discontinuities. For viscous flows, they 
write the Navier Stokes equations as a 

system of first order equations in terms of 
velocity, pressure and vorticity and then a 
full least squares procedure is applied. The 
resulting algebraic equations are 

symmetric positive definite, but the 
number of unknowns are increased ( 
almost doubled for three dimensional 
problems). 

Needless to say, more work is required 
to determine the optimal form of the 

modification terms and the associated 
solution procedures for the general flow 
case. 


Numerical Results 

Three test problems are solved using 
the present formulation. The first viscosity 
method is used for the first two problems. 
Since the Re is relatively low, no artificial 
viscosity in the momentum or the energy 
equation is needed ( e 2=^=0). For the 

modified continuity equation, a numerical 

boundary condition, =0 at the wall is 

enforced at the wall. The third problem is 
solved by the two viscosity methods. In the 
following some preliminary results are 
presented. 

I. Driven Cavity with A diabatic and 

Isothermal Walls 

Incompressible (M =0) and 

compressible flows ( M m = 0,4) are simulated 

for Re=100 with adiabatic walls. The 
pressure contours of the converged 
solutions are plotted in figures (I-a) and 
(I-b) respectively. Next, the temperature at 
the upper and lower walls are fixed and the 
calculations are repeated. The pressure 
contours are plotted in figures (I-c) and (I- 
d) and the corresponding temperature 
contours are shown in figures (I-e) and (I- 
0. 

The effects of compressibility and the 
difference of wall temperatures, are 
clearly depicted in these figures. 

II- Incompressible and Supersonic Viscous 

Flows over a Wall with and Without a 

groove 

A uniform stream of Re=1000 over a 
flat plate is simulated with the same code. 
The pressure and velocity profiles at 
different locations are plotted in figures 
(Il-a) and (Il-b) for the case of M m = 0. For 

supersonic flow with M m = 3 , the pressure, 

velocity and temperature profiles are 
shown in figures (II-c), (II-d) and (Il-e). 
The pressure, temperature and density 
contours are given in figures (Il-f). (Il-g) 
and (Il-h). As expected, an oblique shock 



wave is formed due to the boundary layer 
displacement effect. 

Next, the calculations are repeated for 
a flat plate with a groove. The pressure 
profiles for the incompressible flow case is 
shown in figure (Il-i), while the pressure, 
temperature and density contours of the 
supersonic case are plotted in figures (II- 
j), (Il-k) and (II-l). 

Ill- Inviscid and Viscous flow Simulations 
of Shock Reflection from a flat plate. 

First, an inviscid supersonic flow ( 
M =2.0) with a reflected shock is calculated. 

oo 

The results are in agreement with the 
exact solution. The pressure contours from 
the two artificial viscosity methods are 
plotted in figure (Ill-a) and (Ill-b). For a 
viscous supersonic flow at Re= 296000 and 
M = 2.0, the results are in agreement with 

OO 

those of MacCormack [15] and with 
experimental data. The velocity profiles 
before, within and after the separation 
bubble are plotted in figure (III-c). The 
skin friction distribution is shown in 
figure (III-d). The pressure contours from 
the first and second artificial viscosity 
methods are compared in figure (III-e) and 
(Hl-f). 

Cartesian grids and bilinear elements 
are used for all the problems tested in this 
paper. Applications to transonic flows over 
airfoils using unstructured finite elements 
are reported in a separate paper [16]. 

Conclusions 

A unified approach for a general flow 
simulation is presented. It is shown that all 
the normalized variables used in the 
formulation are always bounded and the 
proper variation of the density due to 
changes in pressure and temperature is 
recovered in the limit of zero Mach 
number. 

Two artificial viscosity methods are 
applied to obtain numerical results for 
some test cases. The governing equations 
are modified by either Laplacian's or 
Poission's equations for the pressure, the 
velocity components and the total 
enthalpy. Acceptable solutions are 


obtained via a standard Galerkin finite 
element procedure, using equal order 
interpolations for all normalized variables. 

The unified approach, offers a 
convenient formulation which allows, 
using the same code, the simulation of 
compressible and incompressible flows 
where the walls are adiabatic, with 
specified temperature distributions or of 
mixed type. In particular, the low Mach 
number approximation as well as the 
Boussinesq approximation (of the 
buoyancy effects) are special cases of the 
present formulations. 

Obviously, it is always possible to have 
more efficient flow simulations for some 
special cases. For example, when the speed 
of sound is finite, explicit schemes can be 
used to integrate the time dependent gas 
equations, in contrast to the 
incompressible flow case where a 
Poission's equation of the pressure has to 

be solved, at each time step, to guarantee 

the conservation of mass during the time 
evolution process. Another example is the 
special case of incompressible ( constant 
density) flow where the energy 
(temperature) equation decouples and the 
solution of the continuity and the 
momentum equations provide the pressure 
and the velocity components. Therefore, 
the use of the present unified approach for 
the above two examples is more costly 
compared to the use of two separate codes 
tailored for the specifics of these, two cases. 
It is still necessary however to have a 
general code for all speeds and all possible 
boundary conditions to handle the cases of 
mixed nature. 
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Fig.(I-a) Pressure contours 
Incompressible Flow 
Grid 36x36 
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Fig.(I-b) Pressure Contours 
Compressible Flow (M =0.4) 

Grid 36x36 
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Fig.(I-c) Pressure Contours 
Incompressible Flow 
Grid 36x36 

Upper Wall Temp.= 1.0 
Lower Wall Temp.= 0.5 
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Fig.(T-d) Pressure Contours 
Compressible Flow (M r =0.4) 
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Upper Wall Temp.- 1.0 
Lower Wall Temp.= 0.5 
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Fig. (I c) Temperature Contour'; 
Incompressible Flow 

Grid 36x36 

Upper Wall Temp.= 1.0 
Lower Wall Temp.= 0.5 
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Fig.(l-f) Temperature Contours 
Compressible Flow (M f =0.4) 
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Upper Wall Temp.= 1.0 
Lower Wall Temp.— 0.5 
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Fig.(U-a) Pressure r-». 
Incompressible Flow 

Grid 26x31 







Fig.(H-e) Temperature Profiles 
Supersonic Flow (^^,=3) 
Grid 26x31 
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Fig.(IH-a) Pressure Contours 
Inviscid Flow (M^=2) 

Grid 34x44 
Viscosity Method #1 



Fig.(III-b) Pressure Contours 
Inviscid Flow (^4^=2) 

Grid 34x44 
Viscosity Method #2 
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Fig.(lII-e) Pressure Contours 
Viscous Flow (M^l) 

Grid 34x44 
Viscosity Method #1 



